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Abstract 

The classical Helmholtz problem is applied for modelling and numerical investi- 
gation of inviscid cusp-ended separated flow around circular cylinder. Two coordi- 
nate systems are used: polar for initial calculations and parabolic as topologically 
most suited for infinite stagnation zone. Scaling by the shape of the unknown free 
line renders the problem to computational domain with fixed boundaries. Difference 
schemes and algorithm for Laplace equation and for Bernoulli integral are devised. 
A separated flow with drag coefficient Cx = like the so called "critical" flow 
is obtained. The pressure distribution on the surface of cylinder and the detach- 
ment point compares quantitatively very well with the predictions of the hodograph 
method. 

1. Introduction 



In 1868 Helmholtz |TT| introduces the notion of a flow consisting of a potential and 
stagnant zones matching at priori unknown free boundaries which are tangential discon- 
tinuities and where the balance of normal stresses (the pressure) holds. Kirchhoff JT^| 



came up with the first solution for the ideal fiow around fiat plate when the detachment 
points were known in advance. Later on in the turning of our century, Levi-Civita WM, 



Villat [T^ etc. developed further the hodograph method and demonstrated its application 
to fiows around curved bodies. Satisfying an additional condition for smooth separation 
(called now Brillouin- Villat condition |T^, |20|) Brodetsky [Q] obtained by the hodo- 
graph method approximate solution for the circular cylinder with a parabolic expanding 
at infinity shape of the stagnation zone. 

In the years 40 of the present century with the computer advent it was already possible 
to calculate such class of fiows direct at the physical plane. The first calculations |T2|, 
gave interesting results. Along with the Brodetsky fiow a radically different Helmholtz 
fiow takes place with decreasing stagnation zone which forms at infinity cusp |^ . Because 
of the limitation of computers the shape of the zone was not conclusive. It appears that 
the method of hodograph can also be applied to obtain such a fiow (see, e.g., ||lOl) but 
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only for the case of circular cylinder. We also found such cusp-ended stagnation zones 
^ 1^ by means of difference scheme and confirmed by integral-method calculations 
|18[| . A new interesting solution for the shape of the wake behind the circular cylinder 
was obtained after further modifying of the difference scheme. Preliminary results of this 
study are represented in [Q. The features of the algorithm and our further investigation 
will be discussed here. 

2. Posing the Problem 

Consider the steady inviscid flow past a circle - the cross section of an infinitely long 
circular cylinder. The direction of the flow coincides with the line 6* = 0, vr of the polar 
coordinates and the leading stagnation point of the flow is situated in the point 6 = n. 
Taking into account the symmetry with respect to the line ^ = 0, vr we consider the flow 
only in upper half plane. 

2.1. Coordinate Systems 

The gist of our approach is to make use of two different coordinate systems: the polar 
one (turning out to be ineffective for the case of infinite stagnation zones extending far 
away from the rear end of body) and the parabolic one the latter being topologically 
more suited for solving Laplace equation outside infinitely long stagnation zones. We 
initiate the calculations in polar coordinates switching to parabolic coordinates after the 
stagnation zone has fairly well developed and has become long enough. 

In terms of the two coordinate systems (cylindrical and parabolic) Laplace equation 
for the stream function -0 reads: 

-(rV^,), + 4^,, = 0, or ^^(^^^ + ^^J = 0. (2.1) 

The undisturbed uniform flow at infinity is given by 

V'l™ ~^^ooSin6', or ^1™,™ ~ ^^^oo • (2.2) 

On the combined surface "body+stagnation zone" hold two conditions. The first 
condition secures that the said boundary is a stream line (say of number "zero") 

ilj{R{e),e) = 0, e e[o,7T] or ij{s{T), t) = o, t e {o,oo) , (2.3) 

where R{6), S{t) are the shape functions of the total boundary in polar or parabolic 
coordinates, respectively. Here and henceforth we use the notation Fi for the portion of 
boundary representing the cylinder and r2 - for the free streamline (Fig.l). 

Let 6* and r* be the magnitudes of the independent coordinates for which the de- 
tachment of flow occurs. As far as we consider only the case when the stagnation zone 
is situated behind the body then the portion of r2 which describes the free line of the 
flow is defined as < < or r > r*, respectively. On r2 the shape function R{0) 
is unknown and it is to be implicitly identified from Bernoulli integral with the pressure 
equal to a constant (say, pc) which is the second condition holding on the free boundary. 
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For the two coordinate systems one gets the following equations for shape functions R{0) 
or S{t): 



or 



r=R(9) 

o<e<e* 



q + 



2 I 2 

(7 + T 



T < T < OO 



where g is a dimensionless pressure. 

At the symmetry line 9 = O^n additional conditions are added 



1 . 



(2.4) 



ct=S{t) 



or -I- = , r = . 

OT 



(2.5) 



and thus (|2.1|) , (|2.2| ), ( p.3|) , ( p^) and ( p.5|) complete b.v.p. for stream-function ip. 



2.2. Scaled Variables 



The above stated boundary value problem is very inconvenient for numerical treatment 
mainly because of two reasons. The first is that the boundary lines are not coordinate lines. 
The second is that the shape function of the stagnation zone must be implicitly identified 
from the additional boundary condition ( |2.4| ). Following we scale the independent 
variable [6 or r) by the shape function R{6) or S{t): 

rR-\e) , 



V 



V 



S{r). 



Such a manipulation renders the original physical domain under consideration into a 
region with fixed boundaries, the latter being coordinate lines. In addition the Bernoulli 
integral becomes an explicit equation for the shape function of the free boundary. Scaling 
the independent variable proved very efficient in numerical treatment of inviscid or viscous 
flows with free boundaries (for details see, e.g., P). 

We treat the two coordinate systems in an uniform way denoting ^ = ^ or ^ = r 
depending on the particular case under consideration. In terms of the new coordinates 
(77, (^), the stream function is a compound function i^{'r],6) = ip{r{r],^),^) or ijj{r],T) = 
'0((t(?7, ^), ^) but in what follows we drop the "tilde" without fear of confusion. The 
Laplace equation takes then the form 



where 



a = r] 



V 



R' . 



(2.6) 



or 
I)- 



1, 



c = 



S' 



R 



a = l + S'\ 

with respective boundary conditions (see 

Thus we define a well posed boundary value problem for -0 provided that functions 
R{d) and S{t) are known. On the other hand in the portion r2 of the boundary (where 
these functions are unknown) they can be evaluated from the Bernoulli integral 
which now becomes an explicit equation for the shape function 



R^ + R 



i2 



dr] 



+ R{9) sine 



rj=l 



1 



or 



1 + S" 



o<e<9* , 



r* < T < 00 



dtp 
dr] 



+ T 



r)=0 



1 



(2.7) 
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3. Forces Exerted on the Body 



The presence of a stagnation zone breaks the symmetry of the integral for the normal 
stresses and hence D'Alembert paradox ceases to exist, i.e. the force exerted from the 
flow upon the body is no more equal to zero. Denote by n the outward normal vector to 
the contour F. Then the force acting upon the contour is given by 



R 



- j^pnds = - j( (g + pc)nds = paU"^ [CJ + Cyj] , (3.1) 

where and Cy are the dimensionless drag coefficient and the lifting force. 

After obvious manipulations we obtain for the drag and lifting-force coefficients the 
following expression (see [§) 



C^ = -2 q [R{e) cos 9 + B!{9) sinO] d9 



or 

C„ = 



where the dimensionless pressure is given by 



i?4 



dr] 



+ R{e) sm9 



ri=l 



Cx 

0. 



q [S{r) + S'{r)r] dr 



(3.2) 



1 2 



or q 



1 + S' 



dip 
dr] 



+ T 



ri=0 



. (3.3) 



4. Difference Scheme and Algorithm 



4-1- Splitting scheme for Laplace equation 

For the purposes of the numerical solution, the transformed domain must be reduced 
to finite one after appropriately choosing the "actual infinities". In the case of polar 
coordinates the domain is infinite with respect to coordinate rj only and it fully enough 
to select sufficiently large number r]^ and to consider the rectangle: [0<^<7r;l<77< r]^] 
(Fig. 2a). In the case of parabolic coordinates an actual infinity is to be specified also 
for the r-coordinate, namely Too and to consider the rectangle: [0 < r < Too; 0<r)< r)^] 
(Fig. 2b). In both directions we employ non- uniform mesh. The first and the last r]- 
lines are displaced (staggered) from the respective domain boundaries on a half of the 
adjacent value of the spacing. Thus on two-point stencils second-order approximation for 
the boundary conditions is achieved (see 0. The non- uniformity of the mesh enables us 
to improve the accuracy near the body and to reduce the number of points at infinity. 

In ^-direction the mesh is not staggered but it is once again non-uniform being very 
dense in the vicinity of the rear stagnation point, i.e. in the vicinity of ^ = which is 
of crucial importance when acknowledging the infinity in cylindrical coordinates. It is 
desirable to have the "actual infinity" in cylindrical coordinates as larger as possible in 
order to prepare the ground for switching to the parabolic coordinates. The connection 
between the r-mesh and ^-mesh is derived on the basis of the connections between the 
two coordinate systems, namely 

Tj = ,jR{9j)cos9j + R{9j) , if 0<9j<7r; Sj = j2R{9j) - t] , 

(4.1) 
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and these relations can be transformed when necessary to calculate Sj , tj from Rj , 9j or 
vice versa. 

Due to the topological differences between the polar and parabolic coordinate systems 
after the transition to parabolic coordinates it is necessary to generate a new r-mesh. 
The new mesh has to be sparse at large distances behind the body where the gradients of 
the flow are small. To this end the knots tj are obtained from ( [4.1|) making use of spline 
interpolation. The new r-mesh is uniform on the rigid body and is changing behind the 
body according to the quadratic rule 



r, = {j-l)h, h='-^, j = l,...,[f] + l, 

r, = exp (j -[f]-l)h\nT^, h = ^, j = [f] + 2, . . . , N + 1 



(4.2) 



where [y] is the last point of the rigid body 

We solve the boundary value problem iteratively by means of the method of splitting 
of operator. Upon introducing fictitious time we render the equation to parabolic type 
and then employ the so-called scheme of stabilising correction [^. On the first half-time 
step we have the following differential equations {At is the time increment) 

(4.3) 

for i = 2, - ■ ■ , M, j = 2, . . . ,N with respective boundary conditions 0] 

The second half-time step consists in solving the following differential equations 



I n+1 I 

1~ f ^' = MaAir^% - Ai(aAir)., (4.4) 

for i = 2, . . . , M, j = 2, . . . , N with respective boundary conditions 

Thus the b.v.p. for the stream function is reduced to consequative systems with sparse 
(tridiagonal) matrices (for detail see e.g.,^]. The main advantage of the economical 
schemes of the splitting type is that on each half-time step we solve one-dimensional 
problems with sparse (tridiagonal) matrices. This can be done by means of the Thomas 
algorithm However, the system for streamfunction ipiv^O cannot be solved by 



plane Thomas algorithm since the condition for numerical stability of the elimination 
is not satisfied for all points of domain. For this reason a modification of the Thomas 
algorithm (in fact Gaussian elimination with pivoting for three-diagonal systems) called 
"non- monotonous progonka" (see |Q) is employed for its solution. 



To calculate afterwards the forces acting upon the body we use the simple formulas for 
numerical integration based on the trapezoidal rule, which are consistent with the overall 
second-order approximation of the scheme. 

4-2. Difference Approximation for the Free Boundary 

The equations (|2.7|) can be resolved for the derivatives R'{6) or S'{t) when the radicals 
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exist, i.e. following conditions are satisfied: 



or 



Qir, 



def 



THt) 



> 1 , nr) 



+ R{e) sine 
dip 



dr] 



(4.5) 



T]=0 



where connection between functions T is determined simply by the formula 



m = l{S{r)-S'r)T{r) 



(4.6) 



The above inequalities are trivially satisfied in the vicinity of the leading-end stagnation 



oo or 



point inasmuch as that for ^ ^ vr (or r ^ 0) one has T — > and hence 

oo. In the present work we use the dynamic condition (|2.4]) in polar coordinates 



T2 

only, so that we present here just the relevant scheme in polar coordinates without going 
into the details for parabolic coordinates. 

Suppose that the set functions S"", T"" are known from the previous global 

iteration, say of number a.Q We check the satisfaction of ( [4. 51 ) beginning from the point 
6 = and continue with increasing 6 . Let j* + 1 be the last point where ([4.5|) is satisfied 
and, respectively j* - the first one where it is not (polar coordinates). The position 9* of 
the detachment point is captured by means of a linear interpolation 



9* 



9* — 9j*+i 



For the shape function Rj of free line is solved the following difference scheme 



Rj + R 



■i-i 




+ 



(4.7) 



for j = j*,... ,2 , whose approximation is 0{gj). Only in the detachment point the 
difference scheme is different, specifying in fact the initial ("inlet") condition, namely 



Rj* 



R{9* 



.R{9*) + Rj, 



\ 



Rp 



+ 



'R{9*) 
T{9*) 



(4. 



where R without a superscript or "hat" stands for the known boundary of rigid body. 
Thus the mere condition for existence of the square root of the Bernoulli integral defines at 
each iteration stage a the new approximation for the position of the detachment point so 
that it 'slides' during the iterations alongside the rigid body. This manner of determining 
of the detachment point we called Christov's algorithm (see ||18|| ). 



^We distinguish here between global and local iteration, the latter referring to the time-stepping of 
the coordinate splitting method. 
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In the end a relaxation is used for the shape-function of the free boundary at each 
global iteration according to the formula: 

= ujRj + (1 - u;)i?J 

where uj is called relaxation parameter. 



4-3. The general Consequence of the Algorithm 

Each global iteration contains two stages. On the first stage, the difference problem for 
Laplace equation is solved iteratively either in polar or in parabolic coordinates (depending 
on the development of the stagnation zone). The internal iterations (time steps with 
respect to the fictitious time in the splitting procedure) are conducted until convergence 
is achieved in the sense that the uniform norm is lesser than £2 = 10~^. Thus the new 
iteration for stream function ibf^'^ is obtained. 

The polar coordinates appear to be instrumental only on the first several (7-10) global 
iterations. When the rearmost cusp point of the stagnation zone reaches 30-50 diameters 
of cylinder (calibers), the current-iteration values of the sought functions are transformed 
to parabolic coordinates and hence the calculations for the stream function continue solely 
in terms of parabolic coordinates. 

The second stage of a global iteration consists in solving the difference problem for 
the free surface in polar coordinates. The transition to and from parabolic coordinates is 
done according to ( [4.1|) and (|4.6| ). Note that there is one-to-one correspondence between 
the points in polar and parabolic coordinates and hence between the respective values of 
the scalar set functions if) and R. 

The criterion for convergence of the global iterations is defined by the convergence of 
the shape function as being the most sensitive part of the algorithm, namely the global 
iterations are terminated when 



max 
3 



Rj 



a+l 



< 10"^ (4.9) 



The obtained solutions for the stream function and the shape function of the boundary 
are the values of the last iteration ipj^^ = i/Jij^^ and Rj = R'j^^ , respectively. Then the 
velocity, pressure, and the forces exerted from the flow upon the body are calculated. 



5. Results and Discussion 



The numerical correctness of scheme (4.3), (4.4) is verified through exhaustive numer- 
ical experiments and through comparison with the known exact solution for the inviscid 
non-separated flow past a circular cylinder 

^ = f/oo(r--)sin6', (5.1) 
r 

where ip is the stream function, Uoo ^ the velocity of the main flow and r and 6 - the 
polar coordinates of a point of the flow. We used different meshes with sizes M x N : 
41 X 68, 41 X 136, 161 x 158, 101 x 201, 101 x 136, etc. Respectively, the actual infinity 
?7oo assumed in the numerical experiments the values 5, 10, 20. The dependence of the 



7 



numerical solution on the time increment At is also investigated and it is shown that the 
approximation of the stationary part of the equations (4.3) and (4.4) does not depend 
on At, i.e. the scheme has the property called by Yanenko ||21[ full approximation. The 
relative differences for ip when At is in the interval [0.001,2] do not exceed 0.5%. The 
numerical experiments show that the optimal values for At is in interval [0.5,1]. For this 
reason the rest of the calculations in the present work are performed with At = 0.5. The 
comparison of the solution ( |5.1| ) to the present numerical results is quantitatively very 
good. The deviations for the different meshes are in order of approximation 0{h^ + g'^) 
and do not exceed 3%. For example in case of mesh 161 x 156 the relative error is about 
0.9%. 

The adequate choice of the "actual infinities" r]^,Too and the spacings hi,gj have a 
profound impact on the accuracy of the difference schemes (4.3) and (4.4). For a given 
"actual infinity" the improvement in the accuracy can be achieved through increasing the 
number of mesh points (decreasing the size of spacing). This makes the use of uniform 
mesh ineffective because in the far-field region the gradients of the flow are small and the 
high resolution is not necessary. That was the reason to employ the non-uniform meshes. 
The "optimal" value for the relaxation parameter turned out to be a; = 0.01. Smaller 
values increased intolerable the computational time while a; > 0.1 could not ensure the 
convergence of the global iterations. Respectively rj^ = 10 is the optimal value for the 
lateral "actual infinity" 

In order to compare calculated results with the prescription of the Levi-Civita method 
in case of so called 'critical' separation angle 6^, = 124.2" (in respect to leading stagnation 
point of the cylinder) it is necessary to summarize that method and deduce corresponding 
relations. Following |jl| the physical plane z is mapped on the unit halfcircle t so, that free 
boundary transforms into the diameter and rigid boundary - into the halfcircumference 
t = e^'^,0 <(T<TT. Then 

..-^'...,,_.,(,_^)*, ,,.) 

where the function Q{t) = 6(t) + iT{t) = Y^k'=o '^2k+it'^^~^^ ■ Hence we obtained the follow- 
ing parametrical equations describing boundary of the cylinder from the leading stagnation 
point to the separation point 9^,: 



Xcyi{s) = Rez = — M J e sin 6(t) sin ct(1 + sin cr)d(T 

.s" (5-3) 
Ucyiis) = Imz = M / e ^"-'^ cos0(t) sincr(l + sin(T)(i(T , 

where f < s < tt , Q{t) = J2k=Q <^2k+i cos{2k + l)a , T{t) = Efclo «2fe+i sin(2A; + l)a , 
and parametrical equations describing the freestreamline from the separation point 6** to 
infinity: 



x{s) = Hez 
y{s) = Imz = ^ / 1 



M r 
M 



- t 



t^ 



cose(t) 



-2 sine (t) 



t3 



--t-- sin0(t)+ --2 cose(t) 



t2 



dt + Xcyiin 



dt + y^yiin) 



(5.4) 
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where — 1 < s < , 9(t) = Y^I^q a2k+it'^''~^^ ■ If the parameters have values M = 
5.71464, ai = 2, ag = .12518, as = .02661 , ay = .00858, ag = .00349, an = .00167, aig = 
.00089, ai5 = .00053, ai7 = .00035, aig = .00024, oai = .00018, aaa = .00016, it corresponds 
to the so called critical separated flow, which detaches at angle 6**. This Helmholtz flow 
has decreasing (concave) stagnation zone with cusp end at infinity (Chaplygin-Kolscher 
flow). 

Further the velocity 

v(z) = iltiie-«f^(*) from where \v(z) \ = — — e'^^*^ , (5.5) 

whence it follows immediately that the pressure on the cylinder is 

p[e) = l-\v\^, (5.6) 

where 6 = arctan is the polar angle. 

In Figs.3-a,b are presented the obtained shapes of the stagnation zone behind the 
cylinder and in the near wake for resolutions 41 x 68, 81 x 136 and 101 x 201 and dif- 
ferent values of relaxation parameter: oj = 0.01; 0.001. Obviously there is an excellent 
comparison between different numerical realizations. On the same figure is added the 
shape of the Chaplygin-Kolscher flow. The latter we calculate by means of parametrical 
equations (|5l3|), ( |5.4| ) using the usual trapezoidal rule. The symbols stand for the results 
taken from the charts of the paper [|17|. It is worth noting the perfect coincidence of the 



computed by us separation angle with both this one, computed in and the 'critical' 
one, prescribed by the hodograph (Levi-Civita's) method. Nevertheless the difference 
between our solution and this in |]T^ is sizable due to the inconclusive character of the 
latter. The logarithmic scale is used in Fig.3-b in order to expand the differences be- 
tween the different difference solutions making them visible in the graph. The shapes 
of the free boundary obtained on the three grids with different resolution are compared 
among themselves very well up to 200 calibers . It is clearly seen up to 70 calibers the 
shapes are practically indistinguishable and up to 160-170 calibers the relative difference 
does not exceed 1-3% respectively. This supports the claim that indeed a solution to the 
Helmholtz problem has been found numerically by means of the developed in the present 
work difference scheme. At the Fig.3-b it is seen the quantitative difference between our 
numerical solution of cusp-ended type and this one prescribed by the hodograph method. 
Indeed there is excellent agreement concerning the positions of detachment point and 
pressure distribution but the hodograph method postulates the asymptotic behaviour of 
the free line also. On the contrary we do not set any condition at infinity. In a sense our 
free boundary has an implicit numerical "closure" of cusp-ended type. 

The calculated dimensionless pressure q is shown in Fig. 4. Here is seen again an 
excellent comparison among the different mesh resolutions. In the stagnation zone it is 
in order of 10~^, which is in very good agreement with the assumption that the unknown 
boundary is defined by the condition q = 0. The amplitude of the minimum of q is 
smaller than 3 the latter being the value for ideal flow without separation. This means 
that the stagnation zone influences on the flow upstream. On the same figure is presented 
the pressure, calculated by means of ( |5.6D which corresponds to the separation angle 
9^ = 124.2°. Apparently obtained here pressure approximates very well this curve. It 



9 



is known the Chaplygin-Kolscher flow has a vanishing drag coefficient [jT], [T^])- In other 
words there exists an inviscid separated ffow submitted to The D'Alembert paradox hke 
nonseparated. Varying the mesh parameters we obtained for the drag coefficient Cx values 
between 2 x 10~^ when resolution is 41 x 68 and 5 x 10~^ when resolution is 101 x 201. 
That is to say our ~ and the error is in order of approximation. In order to confirm 
the above assumption we made the following numerical experiment: in formula (|3.2|) for 
the drag coefficient we replaced our pressure by the pressure obtained from ( |5.6| ). The 
calculated value is C^; = 3 x 10""^. 

6. Concluding Remarks 

An algorithm for numerical solving the classical Helmholtz problem behind a circular 
cylinder is developed. Scaled coordinates are employed rendering the computational do- 
main into a region with fixed boundaries and transforming the Bernoulli integral into an 
explicit equation for the shape function. The crucial feature of the method developed here 
is that the detachment point is not prescribed in advance. Rather it is defined iteratively. 
Difference scheme using coordinate splitting is devised. Exhaustive set of numerical ex- 
periments is run and the optimal values of scheme parameters are defined. Results are 
verified on grids with different resolutions. The drag coefficient of the calculated separated 
flow vanishes like cusp-ended infinite fiow obtained by means of the hodograph method. 
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FIGURE CAPTIONS 

Figure 1: Posing the problem 
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Figure 3: The obtained separation lines for relaxation parameter u = 0.01 and different 
resolutions: 41 x 68; 81 x 136; 101 x 201; hodograph method; 
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Figure 4: The pressure distribution for relaxation parameter u = 0.01 and different 

resolutions: 41 x 68; 81 x 136; 101 x 201; hodograph method; 

o o o o nonseparated inviscid flow. 
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